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DNA cyclization is a powerful technique to gain insight into the nature of DNA bending. The 
worm-like chain model provides a good description of small to moderate bending fluctuations, but 
some experiments on strongly-bent shorter molecules suggest enhanced flexibility over and above that 
expected from the worm-like chain. Here, we use a coarse-grained model of DNA to investigate the 
thermodynamics of DNA cyclization for molecules with less than 210 base pairs. As the molecules get 
shorter we find increasing deviations between our computed equilibrium j-factor and the worm-like 
chain predictions of Shimada and Yamakawa. These deviations are due to sharp kinking, first at 
nicks, and only subsequently in the body of the duplex. At the shortest lengths, substantial fraying 
at the ends of duplex domains is the dominant method of relaxation. We also estimate the dynamic 
j-factor measured in recent FRET experiments. We find that the dynamic j-factor is systematically 
larger than its equilibrium counterpart, with the deviation larger for shorter molecules, because not 
all the stress present in the fully cyclized state is present in the transition state. These observations 
are important for the interpretation of recent experiments, as only kinking within the body of the 
duplex is genuinely indicative of non-worm-like chain behaviour. 


INTRODUCTION 

As the mechanical properties of DNA play an important 
role in its biological capacities, there has been much activ¬ 
ity to accurately characterize these properties, not only 
in the elastic regime of small to moderate fluctuations 
but also for more strongly stressed systems. For example, 
DNA is found to overstretch beyond a salt-dependent 
critical force [1]. Similarly, in response to twist DNA 
forms plectonemes beyond a critical buckling superhelical 
density [2]. Here, we are interested in the response of 
DNA to strong bending. 

The worm-like chain (WLC) model provides a good 
description of small to moderate bending fluctuations in 
DNA [SHS]. However, although there is a consensus that 
for sufficiently strong bending the stress will be localized 
within small regions, often termed “kinks”, much about 
this crossover to non-WLC behaviour remains controver¬ 
sial. For example, a recent review by Vologodskii et al. [7] 
highlighted a number of open questions, including what 
is the free energy cost of kink formation, how does the 
free-energy of a kink depend on bend angle, what is the 
critical curvature that causes the double helix to kink? 

One way to address these questions is with molecular 
simulations of a coarse-grained DNA model, as such sim¬ 
ulations are able to directly probe the relevant factors. 
In this and an accompanying paper [8] we have begun 
this task for the oxDNA model isHm, which is able to 
describe well the thermodynamics of DNA hybridization 
and basic mechanical properties such as the persistence 
length and the torsional modulus [9] and which has been 
productively applied to study a wide variety of biophys¬ 
ical and nanotechnological systems [12]. In the first paper, 
we addressed the thermodynamics of kink formation and 


one particular experimental system, a “molecular vice” 
m, that probed DNA strong bending. Here, we consider 
DNA cyclization, in particular looking at the length of 
DNA molecules at which duplex kinking begins to play a 
role in this process. 


DNA Cyclization 

DNA cyclization is a convenient model system used to 
probe DNA bending. Cyclization experiments were first 
reported in 1966, albeit on 48 500 base pair (bp) A-DNA 
HU HD. In 1981, Shore et al developed a method to 
probe the bending of shorter 126-4361 bp fragments m, 
later noting periodicity in the cyclization efficiency of 
237-254 bp fragments [3|. 

More recently, there has been a particular interest in 
probing the cyclization of sub-persistence length DNA, 
to explore whether this regime is accurately described by 
the WLC model. For example, in 2004 Cloutier & Widom 
(C&W) [T7| challenged the conventional wisdom of WLC 
flexibility established by Shore et al [3 US], claiming much 
greater than predicted cyclization efficiency in 93-95 bp 
DNA fragments. This controversial finding spurred debate 
on the characteristic length at which DNA cyclization 
efficiency deviates from the predictions of the WLC model. 
Despite much experimental la [iMi] and theoretical 
effort [71 I22H35] . a consensus has not yet been established. 

A typical cyclization experiment, as depicted in Fig¬ 
ure [^a), uses a cyclization substrate with complementary 
sticky ends, TVg bases in length, on both ends of a 
base-pair duplex. Cyclization leads to the formation of 
a TVbp-base-pair duplex, where TV^p = TVg + The 
resultant structure is not a closed minicircle - two back- 
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(a) 


kcyc 

^uncyc 


(b) 


kdh 


^undim 



Figure 1. Schematic representations of (a) cyclization where 
kcyc and A;uncyc are the forward and reverse rate constants 
respectively, and (b) dimerization where the rate constants are 
^dim and A^undim- Note that for the dimerization system there 
is only one complementary sticky end per monomer, the other 
end being blunt to allow for only one reaction product, a linear 
dimer. Figures are oxDNA representations for monomers of 
length Abp = 101, including complementary sticky ends of 
length Ns = 10. (Dimer length is 2iVd + iVs). 



bone ‘nicks’ are present at either end of the sticky ends. 
Either the forward rate or equilibrium constant of the 
cyclization reaction is reported. Experiments differ in 
how exactly they probe cyclization: methods based on 
ligation [3l lU [161 , FRET [T9]-[2T] and multimerization 
[36] have been reported. 

In ligase-based experiments, cyclized molecules are 
fixed in the cyclized state by ligation of the two backbone 
nicks. The open and ligated cyclized molecules can then 
be resolved by gel electrophoresis and the concentration 
of different products measured. FRET-based experiments 
can be performed in equilibrium, with the molecules al¬ 
lowed to cyclize and uncyclize indefinitely. Fluorophores 
are attached to both ends of the molecule as FRET report¬ 
ers: a high FRET signal will be reported when the duplex 
ends are in close proximity (cyclized), low FRET when 
apart (open). Although non-WLC behaviour has been 
suggested by the ligase-based experiments of C&W m 
and the FRET-based experiments of Vafabakhsh & Ha 
(V&H) dS], these results have been contested [iiiiET]. 
There is not yet a consensus on their interpretation. 

Cyclization efficiency is traditionally reported in terms 
of a j-factor, first introduced by Jacobson & Stockmayer 
[38] , which is a measure of the effective local concentration 
of duplex ends at zero end-to-end separation. The j-factor 
enables the ring closure probability to be calculated, and 
importantly, may be related to a ratio of equilibrium 
constants: 


K^yc 

i-factor=jeq= (1) 

where are the equilibrium constants for 

cyclization and dimerization, respectively. 

Multimerization of the cyclization substrate yields a 
mixture of linear and circular products [siin]. To avoid 
this complication, a separate dimerization substrate may 
be prepared with the same sequence as the cyclization sub¬ 
strate, but with only one TVg-base complementary sticky 


end per molecule uni (Eigure[^(b)). Following hybrid¬ 
ization, the total length of the system is then Ng 

base pairs, with blunt as opposed to sticky ends. The 
consequences of this choice for and hence jeq are 

discussed in Supplementary Section fSl C[ 

By assuming the contribution from base pairing to Ag^^ 

and A^^ is the same, the WLC model can be used to 
estimate jeq- A common assumption is that the cyclized 
state is fully stacked, with coaxial stacking across the two 
nicks. For this situation, the analytic expression derived 
by Shimada & Yamakawa (SY) (3^, which includes both 
the bending energy cost of bringing the two ends together 
and the twist energy cost of bringing the two helix ends 
into register, is the most appropriate. 

Some of the confusion surrounding claims of cyclization 
efficiency greater than that predicted by the WLC model 
revolves around the use and interpretation of j-factors. 
While the j-factor relation using the ratio of cyclization 
to dimerization equilibrium constants is well established, 
experimental studies usually report the ratio of forward 
rate constants. In the case of ligase-based assays (reviewed 
in reference na): 

.ligase 

• hgase _ '^cyc 

Myn L ligase ’ ^ ' 

Mim 

where and k^^^ are the forward rate constants for 

the formation of the ligated circle and dimer, respectively. 

In the experimental limit where the ligation rate is very 
slow compared to the rate constants for uncyclization 
(^uncyc) and undimerization (A:undim )5 iFe concentrations 
of unligated circles and dimers will reach an equilibrium 
with reactants. If this condition is met, should 

be equivalent to jeq- fa practice, this limit is valid for 
low ligase concentrations ([Ligase] ^ [DNA]) and short 
sticky ends 0 US]. The importance of this condition 
is illustrated by Du et al. [5], who suggested that the 
apparent non-WLC behaviour in the C&W experiments 
was rather due to an insufficiently low ligase concentration 
(reviewed in reference [7]). 

An additional complication with ligase experiments 
relates to the structure of the substrate, specifically the 
DNA surrounding the nick. It is unclear whether the 
ligase will act equally on all nicked duplexes, or only the 
subset that happen to adopt the right configuration at 
the nick, be that coaxially stacked or kinked, to allow 
the ligase to bind. As we show in this work, to alleviate 
stress, cyclized systems are more likely than dimers to 
break coaxial stacking at a nick. The ligation rate of 
hybridized complementary sticky ends may therefore vary 
substantially depending on the system. 

The FRET experiments of V&H [T9| have the advant¬ 
age of directly monitoring the transition between cyclized 
and open states. While it is possible to report thermody¬ 
namics from FRET experiments, V&H report dynamics: 


•FRET ^ ^cyc 


(3) 
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where kcyc and are the forward rate constants for the 
formation of the unligated circle and dimer, respectively. 

Dynamic j-factors extracted from FRET-based exper¬ 
iments must also be interpreted with care. Making a 
comparison between and the SY prediction for 

V&H make a claim of much greater than WLC 

flexibility at TV^p < 100 bp. However, ~ ieq only 

in the limit where A:uncyc ~ ^undim? a cond.ition that V&H, 
as well as another more recent FRET experiment m, 
have shown not to be met. Given that ^uncyc 7^ ^undim^ 
one should not expect ^ jeq at short TV^p. Thus, 

the observed deviation of from is not neces¬ 

sarily an indication of non-WLC behaviour. 

The interpretation of cyclization experiments is not 
straightforward; in particular, the microscopic states, re¬ 
sponsible for the putative non-WLC flexibility, cannot be 
directly observed. Additionally, non-WLC behaviour has 
been reported in a number of other systems, including 
DNA minicircles gn im, a “molecular vice” di, and 
a “strained duplex” [i3Hi5] . As we emphasised in the 
accompanying paper [8], establishing whether these ob¬ 
servations are mutually consistent is a challenging task 
due to the distinct interplay of mechanics, geometry and 
topology inherent in each experimental system. 

Fortunately, simulation can help bridge this gap. Mod¬ 
els that reproduce the basic mechanical, structural and 
thermodynamic properties of DNA can be used to sim¬ 
ulate the relevant systems. The results can be used to 
establish: (i) whether experimental results are truly indic¬ 
ative of non-WLC behaviour; (ii) whether the behaviour 
can be attributed to certain types of structure, such as 
‘kinks’ [TOf48] : (iii) whether different systems present con¬ 
sistent evidence for non-WLC behaviour; and (iv) whether 
the inferred occurrence of any disruptions to duplex struc¬ 
ture is consistent with our current understanding of DNA 
thermodynamics. 

Here we use oxDNA to probe the cyclization equilibrium 
as a function of duplex length. This approach provides 
direct access to microscopic states, enabling a structural 
interpretation of experimental observations. OxDNA is 
particularly well-suited to this task, as it provides a good 
description of both single- and double-stranded DNA, in¬ 
cluding hybridization thermodynamics, persistence length, 
torsional modulus and basic structure, and has previously 
been shown to reproduce a number of stress-induced trans¬ 
itions in DNA [49H52] . 

MATERIALS & METHODS 
oxDNA model 

OxDNA [9UTT] is a nucleotide-level coarse-grained 
model of DNA that has been employed successfully for a 
wide variety of systems m, beginning with the thermody¬ 
namic and structural characterization of DNA nano tweez¬ 
ers [53]. Briefly, the model consists of rigid nucleotides 


with three interaction sites per nucleotide, interacting 
via Watson-Crick base pairing, base stacking, excluded 
volume, and a potential to represent backbone connectiv¬ 
ity. The model is parameterized to reproduce the thermo¬ 
dynamics of duplex melting at high-salt ([Na^]=500 mM), 
where backbone-backbone electrostatic repulsion is short- 
ranged due to counter-ion screening. Two parameteriza- 
tions of the model are available. We use the the average- 
base parameterization pun], in which the strength of 
the base-pairing and stacking interactions are independ¬ 
ent of the identity of the bases, to highlight the basic 
thermodynamics of DNA cyclization. By contrast, the 
sequence-dependent parameterization m has stacking 
and base-pairing interactions that depend on the base 
identity and is used to compare more directly to the V&H 
experiments. 


Simulations 

Simulations of the cyclization equilibrium were per¬ 
formed with a virtual-move Monte Carlo (VMMC) al¬ 
gorithm [54| at 298 K. As the free-energy barrier between 
typical open and cyclized states is large, the transition 
between the two macrostates constitutes a rare-event. 
Umbrella sampling, a technique that allows for the biased 
sampling of states with respect to an order parameter m, 
was employed to sample the barrier crossing in reasonable 
computational time. 

We use a two-dimensional order parameter Q = 
(Qee? Qbp) fo characterize the transition. Qee is a discret¬ 
ized measure of the distance of closest approach between 
the complementary sticky ends. Qbp is the number of 
base pairs formed between complementary sticky ends, 
where 0 < Qbp < ^s- 

To further improve computational efficiency, umbrella 
sampling was windowed to separately sample the open 
and cyclized states of each molecule. For the window 
associated with the open state, the system was restricted 
fo Qhp = 0; for the window associated with the cyclized 

state, the system was restricted to Qee = (fhe value 

corresponding to the shortest distances between sticky 
ends). Simulations were run until convergence to within 
±5 % for each window. 

The sampling windows overlap at Q = Qbp)’ 

results were combined by normalizing each window so 
that the free energies were equal for this value of the 
order parameter. As there is only one well-defined over¬ 
lap between the values of the order parameters for both 
windows, more complex approaches, such as the weighted 
histogram analysis method ISSEZI, were unnecessary. To 
further simplify sampling, we forbade the formation of 
base pairs that are not intended in the design of the 
system (non-native base pairs). Further details of the 
cyclization simulations can be found in Supplementary 
Section [5 1 A[ The simulation of dimerization equilibrium 
is roughly analogous to cyclization, and is elaborated in 
Supplementary Section [Sl B[ To compute the equilibrium 
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constants, we deemed all states with Qbp > 1 to con¬ 
tribute to th e cyclized and dimerized states (details in 
Section |S1 Cl) . 


A complete list of sequences is available in Table SI 


Error bars represent the standard error of the mean from 
5 independent simulations. 


Structural analysis 

For sufficiently strong bending stress, localized struc¬ 
tural disruptions to the DNA double-helical structure are 
expected. We define three such disruptions, namely fray¬ 
ing, bubble formation and kinking, which are elaborated 
in detail in the accompanying paper [8]. Briefiy, both 
fraying and bubble formation involve the breaking of base 
pairs; the difference is in the location along the duplex. 
Fraying involves disruption of base pairing at the duplex 
ends, while bubbles occur in otherwise fully base-paired 
contiguous stretches away from the duplex ends. 

Conceptually, a kink is an area of strong bending local¬ 
ized to a small segment of DNA, and can occur both at 
a nick and within the duplex. When kinking occurs in a 
duplex region, it is nearly always accompanied by bubble 
formation. Similarly, kinking at a nick can be accompan¬ 
ied by fraying. We do not attempt to distinguish different 
types of kinks, as has been done when analysing atom¬ 
istic molecular dynamics simulations [13 Eg. Instead, we 
simply define a kink as present within a duplex region 
when there is a change in orientation of consecutive bases 
of greater than 90° on either strand. For nicked regions, 
only consecutive bases on the intact (unnicked) strand 
are considered. Although this cutoff is somewhat arbit¬ 
rary, since in the current system, kinks arise to localize 
bending stress, they are usually very sharply bent and 
the criterion works well [8]. Sometimes, however, it gives 
rise to false negatives, particularly in the case of kinks at 
nick sites, and when fraying is present. Nonetheless, it is 
satisfactory as an indicator of behaviour for our purposes. 
More details, along with subtleties related to kinking at 
a nick, are discussed in Supplementary Section |S1 D| 


RESULTS 


We simulate a large range of system sizes; some illus¬ 
trative configurations are shown in Figure Quantitat¬ 
ively, we first consider the behaviour of jeq as a function 
of length. For the dimerization system we only com¬ 
puted the equilibrium constant for a few lengths 

(monomers App = 30,67,73,101). As expected, we found 


to be length-independent to within numerical error 
(Supplementary Section S2A, Table S2); therefore, we 
use an average value of = 0.92 ±0.20 x 10^^ in 


our j-factor calculation. In contrast, we found and 

thereby jeq (Eq. 0)> to vary substantially with length. 
In Figure we show jeq^^^ values, calculated from 


our measured values of A^^^ and A^™ using Eq. 
for 81 different lengths in the range A^p = 30 — 2t 
for fixed Ag = 10 using the average-base parameteriz¬ 
ation of oxDNA. These results are compared to 
predictions based on the Shimada & Yamakawa (SY) ex¬ 
pression [39] using previously calculated values for the 
relevant structural and mechanical properties of oxDNA 
[9| . The SY expression is appropriate for the formation of 
a fully stacked “circle” configuration with coaxial stacking 
at both nicks (Figure [^(a)). Note that the comparison 
in Figure is fit-free. The effect of varying Ag as well as 
the role of nicks and mismatches are discussed in Supple¬ 
mentary Section |S2 C 


The behaviour of the SY expression is well understood. 
In the regime of interest, shortening App tends to make 
cyclization less favourable as a result of increased bending 
stress within the duplex. This effect becomes particu¬ 
larly acute for DNA lengths below the persistence length 
(41.82nm/126bp for the curve plotted in Figure]^. On 
top of this systematic behaviour, a periodic oscillation 
is associated with the need to over- or under-twist the 
duplex when the natural twist is not commensurate with 
that required to form a closed circle. The magnitude of 
this oscillation increases at shorter App because a stronger 
twist per base pair is required. 

When considering the behaviour of ieq 
the SY expression, three length-scale dependent regimes 
become apparent: long (A^p > 80), intermediate (A^p « 
45 — 80 bp) and short (A^p < 45). 

In the long length regime (A^p ± 80), oxDNA repro¬ 
duces the periodic oscillations predicted by the SY ex¬ 
pression, and values of jeq coincide at the maxima of 
these oscillations. However, even at the longest lengths 
we consider, the magnitude of the oscillation is smaller 
for oxDNA than predicted by the SY expression. At 
shorter A^p, the magnitude of the oscillation decreases 
for oxDNA, at odds with the SY expression. The differ¬ 
ence stems from the possibility of adopting alternative 
“teardrop” configurations (Figure |^(b)), in which most of 
the twisting stress, and some of the bending stress, can 
be relieved by kinking at one of the nicks, thus breaking 
coaxial stacking. In these configurations there is still the 
constraint of DNA continuity at the nick and this can 
now be more easily satisfied, rather than through over- 
or undertwist, by out-of-plane bending (Figure |^(d)). 

The possibility of adopting this alternative teardrop 
configuration reduces the free-energy penalty for incom¬ 
mensurate values of A^p (i.e. (n ± 1 / 2 ) x pitch length), 
thus suppressing the oscillations. Further, at shorter A^p, 
the bending stress increases more slowly for these kinked 
teardrop configurations than for the coaxially stacked 
circles. Thus, the difference in jeq between the “on-register” 
(coaxially stacked) and “off-register” (kinked at a nick) 
molecules decreases at shorter A^p, rather than increasing 
as predicted by the SY model. Although not taken into 
account for most analyses of cyclization, the possibility 
of the cyclized molecule exhibiting a “teardrop” config¬ 
uration has previously been suggested by Vologodskii et 
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Figure 2. OxDNA representations of different cyclized configurations. Kinks in the duplex, which disrupt stacking and induce 
a 1-3 bp bubble, are indicated with an arrow. All configurations have Ns = 10. (a) A fully stacked “circle”, (b) A “teardrop” 
configuration with a kink at one of the nicks, (c) A = 1 “transition state” configuration, (d) Teardrop configurations 
with Abp ^ {n -\- 1/2) X pitch length can reduce the stress associated with chain continuity at the nick by out-of-plane bending. 
Configurations with a kink in the duplex and either (e) a kink at one of the nicks or (f) kinks at both nicks. For short duplexes 
where is not that much larger than As, the sticky ends can associate either by (g) relatively minor bending of the duplex or 
(h) fraying a few base pairs. 


al [T]. 

Further evidence in support of this analysis is given in 
Figure |^(a), which shows the probability of kinking as 
a function of A^p. At long lengths (A^p > 80), kinking 
does not occur in the duplex regions, but can occur at 
the nick sites. There is clear periodicity in kinking at a 
nick. For example, at A^p = 145 ^ lAx pitch length, the 
probability of kinking at either of the nicks is negligible 
and the system virtually always adopts a coaxially stacked 
circle configuration; however, at the longer A^p = 201 ~ 
19.5 X pitch length, the probability is ^ 40%. As A^p is 
shortened, the probability of kinking at a nick gradually 
increases for these ‘bff-register” lengths; it is not until 
Abp = 114 « 11 X pitch length that the bending stress 
along the duplex is sufficient to cause ^ 10 % of molecules 
to kink at a nick for an “on-register” length. 

In the intermediate-length regime (A^p ~ 45 — 80bp), 
we observe enhanced cyclization efficiency compared to 
the SY prediction. Although continues to de¬ 

crease with Abp, it does so more gradually than the SY 
expression would predict; consequently, is in ex¬ 

cess of the peak envelope of the SY expression. In this 


regime, bending stress in the circular coaxially stacked 
configuration is sufficiently large that kinking at one of 
the nicks occurs for even the on-register systems (Fig¬ 
ure |^(a)). 

Oscillations due to on-register effects also seem to con¬ 
tribute to i°q in the upper end of the intermediate- 
length regime, with shallow maxima occurring at Abp = 
63 and Abp = 74, approximately 6 and 7 times the pitch 
length. However, these sizes no longer correspond to min¬ 
ima in the probability of kinking at a nick (Figure |^(a)) 
so the structural underpinnings of these variations is 
less clear. For shorter lengths, although there are size- 
dependent variations in (e.g. maxima at Abp = 43 

and 49), there is no longer a simple relationship to the 
pitch length, instead reffecting more complex geometric 
compatibilities that allow cyclized states at these lengths 
to be particularly stable compared to nearby lengths. 

In the intermediate regime, most cyclized molecules are 
kinked at one of the two nicks. At shorter Abp, the bend¬ 
ing stress in the duplex region of the teardrop configura¬ 
tions increases, with the highest curvature being localized 
opposite the nick that is kinked. Consequently, it becomes 
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Figure 3. Measured values of for the oxDNA average-base parameterization (black circles) as a function of Abp for 

Ns = 10. For comparison, the Shimada & Yamakawa (SY) WLC prediction [39] (grey solid line) is plotted using values of 

torsional stiffness (4.75 x Jm~^), persistence length (41.82nm) and pitch length (10.36 bp/turn) appropriate to oxDNA 

EllSl]. The dashed line gives the maxima envelope for the SY prediction (grey dashed line). Sequence-variation in 
computed using oxDNA’s sequence-dependent parametrization, is illustrated for sequences considered by V&H (brown triangles), 
including six at A^p = 73 bp, as well as poly-AT and poly-GC. 


increasingly favourable to localize bending stress into a 
kink in the duplex (unnicked) region, with the probability 
of this duplex kinking (Figure |^(e)) increasing from near 
zero at A^p = 81 to near one at A^p = 43 (Figure |^(a)). 
A kink in the duplex will generally be located opposite 
a kink at a nick because this arrangement minimizes the 
residual bending stress in the unkinked portions of the 
duplex by equalizing the lengths of the double-helical 
segments between the two kinks. In addition to a loss of 
stacking, kinking in the duplex typically involves break¬ 
ing 1-2 base pairs [8]. A typical configuration for such 
a kinked duplex state is shown in Figure |^(e). As with 
kinking at the nick, the availability of this configuration 
lowers the free-energy cost of cyclization, raising jeq 
further above the SY prediction. 

At the lower end of the intermediate-length regime, the 
cost of bending without kinking in the duplex region is 
so high that molecules with duplex kinking dominate the 
cyclized state. States with two kinks localize the majority 
of the bending stress at the kinking sites, as shown in 
Figure [^(e), with long duplex sections relatively relaxed. 


Consequently, the free-energy cost of looping is largely 
independent of A^p in this regime, causing to 

level off and reach an approximately constant value. In 
contrast, the SY prediction decreases very rapidly. 

For example, at A^p = 43, ieq^^^ times greater 

than the SY prediction. 

There exists a rich landscape of structures informed 
by slight differences in local geometries as a function of 
Abp. In addition to the canonical two-kink structures, 
containing one kink in the duplex and one kink at a nick 
(Figure|^(e)), we observe several non-trivial arrangements, 
albeit with relatively low probability. For example, the 
configuration in Figure [^(f) contains a kink at both nicks 
as well as in the duplex. 

In the short-length regime (A^p < 42), we observe 
an increase in in stark contrast to the rapidly 

decreasing predicted by the SY expression. Given 

the many WLC assumptions that are violated at this 
length scale, a deviation is unsurprising; however, we did 
not anticipate an increase in 
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Figure 4. (a) Probability of kinking Pkink as a function of 
length in the duplex, at either nick and at both nicks, 

(b) Probability of broken base pairs Pbroken-bp as a function 
of length A/bp for fraying (base-pairing disruption at either 
nick) and bubble formation (base-pairing disruption in the 
duplex region). 


As A/(j is now not that much larger than Ag, kinking 
at both nicks allows the single-stranded sticky ends to 
hybridize without duplex kinking. The system now gener¬ 
ally adopts a conformation of two parallel duplexes, with 
the stress now borne by a mixture of continuous bending 
(Figure [^(g)) and fraying of a few base pairs at the ends 
of the duplexes (Figure [^(h)). Figure [^(a) shows a very 
clear crossover to cyclized states with kinks at both nicks 
and no kinks present in the duplex, occurring abruptly at 
Abp = 40 — 42. At the same time, fraying at the duplex 
ends also increases (Figure [^(b)). In this regime, as Abp 
is shortened, the difference m length of the duplexes de¬ 
creases and the stress in the system tends to drop, leading 
to higher jeq (Figure]^. Overlaid on this overall trend are 
non-trivial geometric effects associated with whether the 
lengths of the two duplexes relative to the pitch length 
is convenient for connecting them (the thickness of the 
double helix is now significant compared to the duplex 
lengths), which leads to non-monotonic behaviour of jeq 
and kinking. Note that the value of Abp at the crossover 
between the short- and intermediate-length regimes is ex¬ 
pected to be very dependent on Ag, occurring at smaller 
Abp for smaller Ag. 


So far we have only reported results using the oxDNA 
average-base parameterization, in which the strength of 
base pairing and stacking interactions are independent 
of base identity. As the free-energy cost of disrupting a 
duplex to form a kink is sequence-dependent, it is import¬ 
ant to consider how sequence might perturb the general 
trends we have elaborated thus far. We therefore studied 
a variety of sequences used in the V&H m experiments 
including six at Abp = 73 using the parameterization of 
oxDNA that includes sequence-dependent thermodynam¬ 
ics (Figure [^. 

As duplex kinking for oxDNA typically involves the 
breaking of base pairs, we observe that this kinking pref¬ 
erentially occurs at A-T base pairs. Base pairs that 
are weaker than average introduce preferred locations for 
kinking; thus, when duplex kinking is relevant, 
is expected to be larger for the sequence-dependent than 
for the average-base parameterization. This is indeed the 
case for Abp = 73. The sequence-induced variation for the 
V&H Abp = 73 sequences is a factor of ^ 4 (GC-content 
11-52%). This compares to a factor of ^ 8 between the 
extrema in GC-content, poly (AT) and poly(GC). 

That we find sequence heterogeneity generally makes 
duplex kinking easier is consistent with our explicit invest¬ 
igation of the free energy of duplex kinking in the accom¬ 
panying paper [8]. Our results imply that the crossover 
to cyclized configurations with duplex kinking occurs at 
slightly longer Abp when sequence-dependence is included 
than for the average-base parameterization. 

As well as equilibrium constants, free-energy profiles as 
a function of the number of bases pairs were computed for 
each system we considered, as these can give more insight 
into the pathway for association. Example profiles for 
cyclization and dimerization at Abp = 101 are illustrated 
in Figure One interesting feature of the profiles is that 
the free-energy gain from hybridizing the complementary 
sticky ends once an initial base pair has formed is less for 


cyclization than for dimerization: AG^ncyc < 

Physically, this indicates that substantial additional bend¬ 
ing stress develops as subsequent base pairs form for the 
cyclization system, reducing the free-energy gain upon 
zippering of the sticky ends relative to dimerization. For 
example, the Qbp = 1 configuration in Figure [^(c) is 
clearly less bent than the Qbp = 1 0 con figuration m Fig¬ 
ure [^(b) (Supplementary Section S2B). The activation 
energies derivable from these profiles will be particularly 
useful in the next section when we consider the dynamic 
j-factor. 


Comparison with experiment 

As noted in the Introduction, the results of ligase 
experiments performed in the low ligase concentration 
limit should, in principle, be comparable to equilibrium 
j-factors (although there are subtleties related to the en¬ 
semble of states actually detected by the ligation enzymes). 
Indeed, for long DNA molecules (Abp much longer than 















Figure 5. Free energy profiles of cyclization (solid) and 
dimerization (dashed) for a = 101 system (Ns = 10, 

= 91). The activation free-energy barriers, AG^ for the 
forward (cyclization, dimerization) and reverse (uncyclization, 
undimerization) reactions are labelled. AGcyc reflects the 
free-energy cost of bending to form the first base pair in a 
cyclization system, whereas reflects the entropic cost of 

bringing two monomers together within the simulation volume. 
The dimerization simulations are for a cubic simulation box of 
dimension 85.18 nm, corresponding to a duplex concentration 
of 2.69 pM. Note that for clarity, we have depicted AGuncyc 
and as the free energy difference between the fully 

base-paired closed state = 10 and the transition state 
Qbp = 1; however, in practice, frayed states Q^p = [2,9] do 
contribute to the closed state. While the distinction does not 
significantly impact our results (< 1/4 /cp ^): we do include the 
contribution of frayed states in both AGuncyc cind 


the persistence length) and low ligase concentrations, 
there has been consistent agreement between experiment 
and the WLC model laiiiiisiiiTiiiniiii]. However, Clou¬ 
tier & Widom (C&W) [17] reported results for TVpp shorter 
than the persistence length (N\^p = 93,94,95,105,116), 
showing an apparent deviation from WLC behaviour, 
with (Eq. 1^) enhanced over the SY WLC predic¬ 
tion [39] by a factor of 10^ — 10^. The differences 

between the maxima in and in this size 

range are much smaller (Figure]^ and so cannot account 
for this discrepancy. 

In contrast, Du et al [5] found no deviation from WLC 
behaviour for Ypp = 105 — 130 . Furthermore, they presen¬ 
ted evidence suggesting that the C&W experiments used 
too high a ligase concentration to enable to be com¬ 


pared with The results of Du et al. are in good 

agreement with the SY WLC expression, albeit with some¬ 
what different materials parameters (torsional stiffness, 
persistence length and pitch length) than for oxDNA, ow¬ 
ing in p art to different buffer conditions (Supplementary 
Section S2D, Figure S4(a)). 


FRET measurements on DNA cyclization, as pioneered 
by Vafabakhsh & Ha (V&H) [T9|, provide a more direct 
measure of cyclization because A;cyc, ^uncyc arid 
obtainable, although V&H mostly report feyc- V&H 
claim enhanced cyclization at App < 100, based on a 

comparison between their (Eq. (j^) and the SY 

WLC expression However, as noted earlier, this 

is only a fair comparison if ^uncyc = ^undim^ which, as 
^undim is expected to be length independent, also implies 
that A:uncyc should be independent of A^p. However, since 
both V&H [19] and more recent FRET measurements [21] 
suggest that A;uncyc increases with A^p, this condition is 
not met. 

Although we do not directly simulate the dynamics 
of cyclization, we can estimate the relative rates of pro¬ 
cesses using free-energy profiles such as those in Figure 
and unimolecular rate theory for activated processes [59] . 
In agreement with experimental investigations, previous 
work on oxDNA has shown that duplex formation has an 
effective “transition state” involving a very small number 
of base pairs m- We therefore make the assumption that 
uncyclization and undimerization rates are given by: 


^uncyc — ^ 


^ndim ^ 


-AG, 


t 

uncyc 


ksT 


-AG, 


ksT 


. \ 

undim | 


(4) 

(5) 


where A is a constant for DNA melting, and AG^^cyc 

AG^ndim defined in Figure]^ The physical content 
of this assumption is that an increased favourability of 
base-pair formation is manifested in slower unbinding 
rates. This is important because bending stress in cyclized 
systems reduces AG^j^^^y^, (Figure [^. The rate constants 
will be equal (A^uncyc = ^ndim) i^ very long length 
limit (Abp much longer than the persistence length). 

AG^yf. and AG^j^^^y^, are plotted as a function of A^p 
in Figure In particular there is a general decrease in 
AG^ncyc shorter lengths, suggesting A^uncyc increases 
with shorter A^p, in agreement with experimental results 
[iai2T]. Additionally, any torsional stress in the cyclized 
state is relieved as base pairs are disrupted, leading to 
the oscillations in AG^^cyc lengths, with the 

minima occurring at the more torsionally stressed off- 
register lengths. 

To compare with V&H’s results for we note 

that 


Jdyn = 


k h 

^cyc ^uncyc 


^dim ^ndim 


( 6 ) 


Therefore, using Eq. rt5|, our approximation for j^y^^ 


oxDNA 
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iVbp 

Figure 6. Free-energy barriers for (a) cyclization AGcyc 
(b) uncyclization AG^ncyc (black) and undimerization (grey, 
average highlighted with dashed line), all computed using the 
oxDNA average-base parameterization. 


is 


•oxDNA 

•^dyn 


eXT) ( ^ ^^ncyc A 

^eq exp \ j 

TZTTd V’ 

^dim / ^^undim \ 

Aeq exp I I 


(7) 


at each App. 

We are now in a position to estimate from 

oxDNA’s equilibrium constants and activation free-energy 
barriers. We expect and to be length- 

independent as the excluded volume of the duplex far from 
the complementary single-stranded sticky ends is likely to 
have little effect on the dimerization process. Indeed, this 
appears to be the case to within less than 0.5 k^T in AG 
(Figure [^b), Supplementary Table ^2 ). Our resulting 
values are plotted in FigureTf, where they are 
compared to our values, and V&H’s 

We observe that oxDNA’s values lie substan¬ 

tially above ieq except in the limit of long App for on- 
register lengths. Only if all the stress-induced destabiliza¬ 
tion of the fully cyclized state is exhibited in the forward 
rate would = jeq- Although the activation free-energy 
for cyclization increases as App decreases (Figure l^a)), 
the full bending and twisting stress present in the cymized 
state is not yet present in the transition state where the 
two complementary sticky ends have formed their first 
base pair, e.g. Figure j^c), and so Iii- 

stead, stress which is not present in the transition state 


leads to a decrease in the free-energy barrier for uncycliz¬ 
ation (Figure [^b)) and accelerates uncyclization relative 
to undimerization. Alternatively, one can consider that 
the stress in the cyclized state subjects the duplex formed 
between the sticky ends to a shear force m which is 
well known to lead to a more rapid rupture of a duplex 

[SI ED. 

The observed behaviour is similar to that seen by V&H; 
indeed, and agree remarkably well in the 

range App = 90 — 105, with reasonable agreement ex¬ 
tending to A^p = 70. Our results suggest that variation 
of uncyclization rates with A^p may be an important 

TjVp IT'X' 

contribution to apparent non-WLC behaviour in 

To begin to understand the difference between V&H’s 
dynamic j-factor and a number of authors have 

tried to account for the role of the single-stranded tails 
in cyclization by incorporating a “capture radius” into 
the WLC j-factor calculation, as an approximation for 
how close the duplex ends must be in order for the sticky 
ends to hybridize m na mi EH- Indeed, this approach 
leads to a significantly enhanced j-factor when a value 
that is taken to be roughly appropriate for the 10-base 
sticky ends of the V&H experiment, namely 5 nm, is used 


naisz]. 

Assuming this approach is an attempt to capture the 
weaker constraints at the transition state, and hence to 
estimate the activation free-energy barrier relevant to j^y^ 
(note jqyn and jeq are often not clearly differentiated in 
discussions of V&H’s results), we can use our oxDNA 
results to test the reasonableness of this approximation 
by measuring the separation of the duplex ends in the 
transition state ensemble. At long lengths our measured 
capture radius is approximately constant with a value 
just under 4 nm, but increases at shorter lengths because 
some of the stress at the transition state is partitioned into 
stretching the single-stranded tails, reaching a maximum 
of about 7 nm at about A^p ^ 40 (Supplementary Fig¬ 
ure S5). Thus, although the 5nm value used previously 
[mi? is not unreasonable, unsurprisingly this approach 
does not capture the full complexity of the transition 
state to cyclization, and nor does it account for the strain 
present in the sticky ends for them to achieve contact. 

Therefore, although explaining a dynamic j-factor in 
terms of a capture radius is physically well-motivated, 
it does not provide a full understanding. Moreover, as 
its value is not tightly constrained, and certain contri¬ 
butions to stability are neglected, it is difficult to judge 
whether WLC behaviour is violated or not by comparing 
a calculated j-factor curve to an experimental Jdyn- 


We note that in contrast to varies 

comparatively smoothly with A^p, with only a weak peri¬ 
odicity on the length scale of the pitch length. The shal¬ 
low maxima at larger lengths occur at A^j = (n + 1 / 2 ) x 
pitch length, because the two sticky ends are then on 
the same side of a torsionally unstressed duplex. Con¬ 
sistent with this smooth variation, ACjy^, also varies 
relatively smoothly with A^p. The strong periodicity in 










10 



Avg-base 
Seq-dep 
^ poly(AT) 
poly(GC) 

^ FRET 
•^dyn 

.^•oxDNA Avg-base 
_I_I_I_I_ I_ 

30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 

^bp 


Figure 7. OxDNA dynamic j-factor (black circle) compared to the FRET experiments of V&H (red squares) 

m- Results for the oxDNA sequence-dependent parameterization using V&H 14 variable sequences, in addition to 6 

sequences at A^p = 73 to highlight the role of sequence-variation. For reference, and are also shown (grey). 


comes from the free-energy gain when zipper- 

ing up the complementary sticky ends {i.e. AG^j^cyc)^ 
which is greater when A^p is an integer multiple of the 
pitch length, allowing the formation of a relatively relaxed 
coaxially stacked circle. V&H suggest that their 
data at A^p = 93 — 106 displays a strong oscillation 
with a period of about one pitch length (Figure [^. In 
agreement with Vologodskii et a/., [7] we find no physical 
mechanism for such a strong oscillation and would suggest 
that the experimental evidence for this oscillation is not 
compelling. 

At the shortest lengths investigated, still lies 

above Some of this difference might be due 

to our use of the oxDNA average-base parameterization. 
Sequence does play a role in DNA flexibility, as noted in 
the accompanying paper [8], but the impact of sequence 
variation is non-trivial. In particular, because kinks tend 
to localize to AT base pairs, kinking in the duplex is easier 
for a sequence with 50 % GC-content than reported for 
our average-base parametrization. 

Using the oxDNA parameterization with sequence- 
dependent thermodynamics we found to increase 
compared to the values for the average-base parameter¬ 
ization. Consistent with our results for we find 

a sequence-induced variation in ^ factor of 

^ 4 for the six V&H A^p = 73 sequences (GC-content 


11-52%) [T9] compared to a factor of ^ 8 between the 
extrema in GC-content, poly (AT) and poly(GC). This 
compares to a factor of ^ 60 for experimental looping 
rates, a discrepancy which may be explained by V&H’s 
use of poly(A) tracts, a sequence-motif well-known to 
introduce intrinsic curvature in duplex DNA [62|. As 
oxDNA’s sequence-dependent parameterization is based 
on the nearest-neighbour thermodynamics of SantaLucia 
et al [631164] . alternative structural motifs such as poly (A) 
tracts are outside the scope of the model. OxDNA also 
does not reproduce sequence-dependent structural (e.g. 
the difference in size between purine and pyrimidine) or 
mechanical (e.g. flexibility) properties. 

The origin of the remaining discrepancy between 

and Jdy?^^ shortest A^p is not yet clear, but may 

indicate enhanced flexibility for V&H versus oxDNA due 
to WLC {i.e. lower persistence length) or non-WLC be¬ 
haviour {i.e. kinking within the duplex at slightly longer 
Abp). The oxDNA persistence length of 41.82 nm is within 
the range of experimental observations at [Na^]=500 mM, 
but smaller values at [Na^]=750mM are not implausible 
[65l l66] . It is also possible that oxDNA slightly underes¬ 
timates the prevalence of kinking within duplex regions 
EH; an onset of kinking at slightly lower stress (longer 
Abp) would make cyclization at shorter lengths more fa¬ 
vourable. Finally, it is worth noting that the presence of 
fiuorophores may cause perturbations in the V&H experi- 














11 


merits. 

V&H do report i^eq^ for some systems; oxDNA results 
are in good agreement for those lengths (Supplementary 
Figure [S^. 

CONCLUSION 

Cyclization is a system dependent manifestation of the 
general thermodynamics of strong DNA bending, elab¬ 
orated in the accompanying paper [5]. The remarkable 
range of behaviour in cyclized systems is explicable by 
the interplay between three specific deformation modes 
of stressed duplexes: continuous bending, kinking and 
fraying. 

OxDNA reveals that each of these modes is present at 
a characteristic length-scale with respect to cyclization: 
continuous bending at long lengths (A^p > 80), duplex 
kinking at intermediate lengths (A^p ^ 45 — 80 bp) and 
fraying at short lengths (A^p < 45bp). In addition, as 
Abp is shortened, there is an increase in kinking at the 
two nicks that remain after the hybridization of the sticky 
ends. At longer lengths, kinking at a nick is only observed 
for “off-register” molecules that cannot form torsionally 
relaxed coaxially stacked circles. The ability of said nicks 
to relax bending, as well as torsional stress, means that 
they become increasingly prevalent for shorter A^p. At 
the shortest lengths, kinking at both nicks is dominant. 

We use oxDNA to probe the reported observation of 
non-WLC behaviour in FRET-based cyclization experi¬ 
ments m- In agreement with experiment, we observe 
that for shorter values of A^p, the apparent j-factor lies 
substantially above the predictions of the Shimada & Ya- 
makawa (SY) [39] WLC model. We also observe that 
the periodic oscillations predicted by the SY model are 
suppressed. This behaviour arises from two conceptually 
distinct phenomena. 

Firstly, highly stressed cyclized systems can adopt con¬ 
figurations that relax stress more effectively than through 
continuous bending, thereby reducing the overall free- 
energy cost of cyclization relative to a direct estimate 
based on a simple WLC-based model. At various val¬ 
ues of Abp, oxDNA identifies kinking at nicks, kinking 
within the duplex region and fraying of base pairs as key 
relaxation modes. 

Secondly, oxDNA suggests that not all of the reduction 
in relative to is due to stress manifest in the 

cyclization rate; uncyclization rates are also substantially 
increased relative to undimerization rates. The result is 
that dynamic j-factors based on the ratio of cyclization 
and dimerization rates lie even further above the SY 
prediction than their equilibrium j-factor counterpart. 

Of the above effects, only kinking within the duplex 
can reasonably be described as truly non-WLC behaviour. 
WLC and related statistical models do not predict ab¬ 
solute rates directly. Kinking at nicks and fraying can 
only occur when the DNA backbone is discontinuous, and 
any resultant effects are unrelated to whether WLC mod¬ 


els accurately describes the body of the DNA duplex. In 
oxDNA (with the average-base parameterization), kinking 
in the duplex only has a substantial effect for A^p < 70 bp. 
Note, if sequence dependence is taken into account, then 
kinking may occur for slightly longer lengths. 

Our results suggest that much of the apparent “extreme 
bendability” reported by Vafabakhsh and Ha m can 
be attributed to factors that are not strictly speaking 
evidence of non-WLC behaviour. We cannot account 
for deviations at their very shortest lengths (A^p < 70). 
In oxDNA, kinking within the duplex region is present, 
but not completely dominant. It is possible that oxDNA 
slightly overestimates the difficulty of kinking within a 
duplex - if this is the case, the data for the very smallest 
values of A^p studied by V&H may be indicative of duplex 
flexibility over and above that predicted by the WLC, 
although we note that to give a substantial effect on 
kinking must not only be present, but must dominate the 
ensemble. Kinking dominating the ensemble at A^p 70 
is inconsistent with oxDNA predictions, and available 
experimental evidence for a “molecular vice” [8| and DNA 
minicircles M- 

Various authors have incorporated a “capture radius” 
into the WLC j-factor calculation to capture phenomen¬ 
ologically some of the effects listed above [Illliia]. In¬ 
deed, this approach leads to an enhanced j-factor, but 
the choice of capture radius is somewhat arbitrary and 
imprecise, making it difficult to assess whether non-WLC 
behaviour is present. Nonetheless, we find that the 5nm 
capture radius that has been used in the interpretation 
of V&H’s experiments is not unreasonable. 

OxDNA is only a model, and good correspondence with 
experimental results should not be over-interpreted. Non¬ 
etheless, the relaxation mechanisms identified are clearly 
physically plausible. For example, enhanced uncyclization 
rates have previously been noted in the literature naiai] 
(and also for transcription-factor mediated looping |68|L 
It is clear that much of the apparent discrepancy between 
the data of V&H and the predictions of WLC-based mod¬ 
els are due to effects that are not true violations of the 
WLC model of duplex DNA flexibility. Our study also 
helps to reconcile the results of V&H with previous ligase- 
based assays which saw no evidence of enhanced flexibility 
at Abp ^ 100 [5], and previous studies of minicircles which 
detected no evidence of duplex disruption at these length 
scales m- 

To explore whether the shortest lengths studied by 
V&H do show evidence of kinking and enhanced flexibil¬ 
ity, we would propose experiments of shorter sequences 
and systematic collection of both dynamic and equilib¬ 
rium data. The latter is extremely important; statistical 
WLC models make equilibrium predictions, and so the 
breakdown of a WLC description can only be confirmed 
with equilibrium data. Indeed, elucidating the subtleties 
of dynamic and quasi-dynamic (C&W) j-factors is one of 
the key issues addressed in this work. 
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Supplementary material for “Coarse-grained modelling of strong DNA bending II: 

Cyclization” 
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SI. SIMULATION METHODS 

A. Cyclization simulations 

Cyclization simulations were performed in three phases: 
exploratory, equilibration and production. In all cases, 
we use a virtual-move Monte Carlo (VMMC) algorithm 
m in combination with umbrella sampling |S2) . 

In the exploratory phase, we iteratively adjusted the 
umbrella sampling bias for the windows associated with 
the open and cyclized states, yielding a flat population 
of states for both windows. We use discrete potentials 
for both dimensions of the order parameter (see main 
text). For Qee, we use variable width distance increments: 
0-1.7036 nm, 3.4072, 5.1108, 8.518, 12.777, 17.036, 21.295, 
25.554, 34.072, 42.59, 51.108 and > 51.108 nm. For Qpp, 
the increment is fixed, and is simply the number of base 
pairs formed (0 < Qpp < Ng). As the weighting itera¬ 
tions were semi-automated, simulation times varied in 
the range 10^-10^ VMMC steps per particle. Simulations 
were performed with one molecule in a cubic box of di¬ 
mension 170.36 nm, which corresponds to a unimolecular 
concentration of 336 nM. To further simplify sampling, we 
forbid the formation of base pairs that are not intended 
in the design of the system (non-native base pairs). 

In the equilibration phase, we equilibrated the system 
for 10^ VMMC steps per particle using the aforemen¬ 
tioned umbrella weights. Adequate sampling (number 
of transitions in Qee and Qbp) and decorrelation (via a 
block average decorrelation method) of potential energy, 
bubble size, fraying and structural kinking were checked. 
For the production phase, each simulation (five independ¬ 
ent simulations per measurement) was initialized with 
a statistically independent starting configuration and a 
unique random seed. This is accomplished by randomly 
drawing starting configurations from the equilibration 
phase, one for each production trial, ignoring the first 
2 X 10^ VMMC steps per particle of equilibration. The 
production trial run time is 10^ VMMC steps per particle. 
For reference, the characteristic decorrelation time for the 
potential energy is ^ 10^ VMMC steps per particle, while 
the decorrelation time for kinking is r\j 10^ VMMC steps 
per particle. 

The “seed moves” used to build clusters in the VMMC 
algorithm were: 

• Rotation of a nucleotide about its backbone site, 
with an axis chosen uniformly on the unit sphere, 
and with an angle drawn from a normal distribution 
with a mean of zero and a standard deviation of 
0.10 radians. 

• Translation of a nucleotide, where the displacement 
along each Cartesian axis is drawn from a normal 


distribution with a mean of zero and a standard 
deviation of 0.085 18 nm. 


B. Dimerization simulations 

Dimerization simulations follow a very similar proced¬ 
ure as cyclization simulations (Section [Sl A[ ), but adapted 
for a bimolecular system. With the exception of the start¬ 
ing configuration (bimolecular), exploratory, equilibration 
and production procedures are identical. As discussed 
in the main text, each molecule has one complementary 
sticky end and one blunt end. This precludes the forma¬ 
tion of multimers and circular dimers, allowing only the 
formation of linear dimers. Simulations were performed 
with these two complementary, non-palindromic, sticky- 
ended duplexes, in a cubic box of dimension 170.36 nm, 
which corresponds to a dimer concentration of 336 nM. 
To estimate the impact of excluded volume effects, dimer¬ 
ization systems were also performed in a box of half the 
size and therefore 8 times the concentration (85.18 nm 
per side, which corresponds to a dimer concentration of 
2.69 pM). As for the cyclization simulations, non-native 
base pairing was disallowed. 

Analogously to the cyclization simulations, the dimer¬ 
ization simulations are windowed using the order para¬ 
meters Qee and Qbp, respectively the distance of closest 
approach and the number of base pairs formed between 
complementary sticky ends. The window associated with 
the dimerized state is given by Qbp > Ibp (Qee = 
the small separation between base sites of the base paired 
nucleotides), while the window associated with the undi¬ 
merized state is given by Qbp = 0 bp. 

C. Computation of equilibrium constants 

The cyclization reaction is unimolecular (Figure [^(a)): 

(51) 

^uncyc 

where A and B are the open and cyclized states respect¬ 
ively, kcyc is the forward rate constant (cyclization) and 
^uncyc is the reverse rate constant (uncyclization). The 
equilibrium constant for cyclization, estim¬ 

ated directly from simulations of a unimolecular, isolated 
system as 

(52) 

Here Pa Pp are the probabilities with which un- 
cyclized and cyclized systems are observed in simulation 

(Fa + = l). 
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The dimerization reaction is a bimolecular association 
of distinct molecules: 


A + B 


^dim ^ 
- 

^ndim 


AB. 


(S3) 


In this case, the bimolecular equilibrium constant can be 
inferred from a simulation of a single pair of dimerizing 
monomers via 


^dim _ ^AB 

"^’a/b[Ao]’ 


(S4) 


in which [Aq] is the total concentration of strand type A 
in simulation, and Pa/B and Pab l^he probabilities 
with which monomers and dimers are observed in simu¬ 
lation. This expression follows from Eq. 7 of Ref. |S4| . 
The resultant obeys the standard relation for bulk 

systems in equilibrium. 


;^dim _ [AB] 

^<1 [A][B]- 


(S5) 


Recall from the main text that we use a definition 
for the equilibrium j-factor jeq that is consistent with 
Vafabakhsh and Ha (V&H) 


Jeq — 


AeT 


^dim ‘ 


(S6) 


This is a measure of the effective concentration of one 
sticky end in the vicinity of the other when the dimer¬ 
ization reaction involves distinct monomers, forming a 
heterodimer, with only one sticky end per monomer. Re¬ 
searchers often estimate using the dimerization of 

identical monomers, forming a homodimer, with palin¬ 
dromic sticky ends. In this case, given the same underly¬ 
ing interaction strength between sticky ends (and hence 
the same A^”^ would be increased by a factor of 

two due to combinatorial (4 times the number of possible 
dimers) and symmetry effects (a factor of 1/2 for a ho¬ 
modimer rather than a heterodimer) |S41IS5| . As a result, 
with this approach the equilibrium j-factor is estimated 
as 

K^yc 

M = 2^- (S7) 


In the case of identical monomers with two non- 
palindromic sticky ends, there is no pre-factor in the 
j-factor estimate, as the combinatorial (twice the number 
of possible dimers) and symmetry effects (a factor of 1/2 
for a homodimer) cancel. Overall, under the approxima¬ 
tion of ideal behaviour of separate complexes, these three 
approaches are equivalent. 


D. Structural criterion for kink detection 

While it is often visually straightforward to identify 
kinks, automating their detection is not without difficulty. 


We have developed both energetic and structural criteria 
for identifying kinks, which are based on disruption of 
stacking interactions and changes in base orientation, 
respectively. In this study, we rely on the structural 
criterion; we have discussed the differences between the 
two criteria in detail in the supplementary material of the 
accompanying paper |S6| . 

Structurally, we define a kink using the relative orient¬ 
ation of consecutive nucleotides. In oxDNA, the orient¬ 
ation of each nucleotide is unequivocally determined by 
two orthogonal unit vectors, the base-backbone vector 
and the base-normal vector. The base-backbone vector 
connects the backbone and base interaction sites of each 
nucleotide. The model is designed such that, in a relaxed 
duplex, the base-normal vector at index z, and at the 
consecutive index i + 1, a^^i, are approximately parallel; 
that is, • a^+i ^ 1. A kink is defined to be present 
if • a^^i < 0, a condition implying a more than 90° 
change in orientation of consecutive nucleotides along one 
strand. Naturally, for duplex regions, we consider kinking 
along either strand; although, usually if a kink is present 
in the duplex, then both strands are kinked. To limit 
false positives due to fraying, we do not include the first 
and last 3 pairs of nucleotides in the duplex region in our 
analysis. 

In the fully cyclized configuration there are two “nicks” 
where the two ends of each strand meet. Kinks in nicked 
regions are treated slightly differently, as our criterion is 
only well defined along an intact strand. A kink at a nick 
is detected on the opposite strand. As kinks at a nick may 
diffuse slightly, we define a region of 3 base pairs on either 
side of the nick on the intact strand. If the intact strand 
is kinked in this 6 base pair region, then the molecule is 
considered kinked at that nick. 

Since distinguishing fraying from kinking at low is 
problematic, we compute kinking only for the most prob¬ 
able values, e.g. = 8 — 10 for the Ng = 10 molecules 
depicted in Figure [4 The lower configurations appear 
in our simulations only because of biased sampling and 
have a negligible contribution to the equilibrium kinking 
probability. 


S2. ADDITIONAL RESULTS 


A. Dimerization equilibrium 

The equilibrium constant for dimerization A^™ is 
length-independent to within two times the standard error 
of the mean for the oxDNA average-base parameteriza¬ 
tion (Table [^. Varying concentration suggests that the 
role of excluded volume effects is small between 336 nM 
and 2.69 pM. 

For the sequence-dependent parameterization, there 
is some variation in A^”^, but this does not represent 
length-dependence per-se, rather it is most likely due to 
sequence variation in the bases at the interface between 
the duplexes and the sticky ends. 
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A/bp Note 
Ns = 10 

30 ^ 

67 

69 

71 

73 

74 
76 

93 

94 
97 
99 
101 
103 

105 

106 
207 


Cyclization sequence 


CAG AAT CCG TGC TAG ACC TCC ACC 
CAG AAT CCG TGC TAG TAG CTC AAT 
CGA TCA TGG AGT TAT ATC TGA GGG 
CAG AAT CCG TGC TAG TAG CTC AAT 
CAG AAT CCG TGC TAG TAG CTC AAT 
CAG AAT CCG TGC TAG TAG CTC AAT 
CAG AAT CCG TGC TAG TAG CTC AAT 
CAG AAT CCG TGC TAG TAG CTC AAT 
CAG AAT CCG TGC TAG TAG CTC AAT 
CAG AAT CCG TGC TAG TAG CTC AAT 
CAG AAT CCG TGC TAG TAG CTC AAT 
CAG AAT CCG TGC TAG TAG CTC AAT 
CAG AAT CCG TGC TAG TAG CTC AAT 
CAG AAT CCG TGC TAG TAG CTC AAT 
CAG AAT CCG TGC TAG TAG CTC AAT 
CAG AAT CCG TGC TAG TAG CTC AAT 

CAG AAT CCG TGC TAG TAG CTC AAT 
TGG GAA AAG AAC GTA GAT GTT AAT 


GTT TCA 

ATA GAC TCC CTT TGA CCT GAC TAT CCT 

AAA CTG GAC TGA TAG GAG TGG AGG TGG 

ATA GAC TCC CTT TGA CCC ATG ACT ATC 

ATA GAC TCC CTT CTA ATT GAC TGA CTA 

ATA GAC TCC CTT CTA ATT GAC CAT GAC 

ATA GAC TCC CTT TAA GTT GAC CCA TGA 

ATA GAC TCC CTT CTA AGT TGA CCT CAT 

ATA GAC TCC CTT AAT ACT TCT CCT ATG 

ATA GAC TCC CTT TAA TAC TTC TCC TAT 

ATA GAC TCC CTA TGT TAA TAC TTC TCC 

ATA GAC TCC CTA GAT GTT AAT ACT TCT 

ATA GAC TCC CTG TAG ATG TTA ATA CTT 

ATA GAC TCC CTA CGT AGA TGT TAA TAC 

ATA GAC TCC CTG AAC GTA GAT GTT AAT 

ATA GAC TCC CTA GAA CGT AGA TGT TAA 

ATA GAC TCC CTA TCA GTA CGA AGC TGG 

ACT TCT CCT ATG ACT TCT AAT TGA CCC 


CAC CTC CAC CGT TTC A 
CAA AGT GTC TTA GGC A 
CTC ACC TCC ACC GTT TCA 
TCC TCA CCT CCA CCG TTT CA 
TAT CCT CAC CTC CAC CGT TTC 
CTA TCC TCA CCT CCA CCG TTT 
GAC TAT CCT CAC CTC CAC CGT 
ACT TCT AAT TGA CCC ATG ACT 
GAC TTC TAA TTG ACC CAT GAC 
TAT GAC TTC TAA TTG ACC CAT 
CCT ATG ACT TCT AAT TGA CCC 
CTC CTA TGA CTT CTA ATT GAC 
TTC TCC TAT GAC TTC TAA TTG 
ACT TCT CCT ATG ACT TCT AAT 
TAC TTC TCC TAT GAC TTC TAA 

GCT ATA CCG TTC TTA TTG TCC 
ATG ACT ATC CTC ACC TCC ACC 


A 

CA 

TTC A 

ATC CTC ACC TCC ACC GTT TCA 

TAT CCT CAC CTC CAC CGT TTC A 

GAC TAT CCT CAC CTC CAC CGT TTC A 

ATG ACT ATC CTC ACC TCC ACC GTT TCA 

CCA TGA CTA TCC TCA CCT CCA CCG TTT 

ACC CAT GAC TAT CCT CAC CTC CAC CGT 

TGA CCC ATG ACT ATC CTC ACC TCC ACC 

TTG ACC CAT GAC TAT CCT CAC CTC CAC 

TTA ATA CCA CCG ACG AGT TGT ACG CCC 

GTT TCA 


CA 

TTC A 
GTT TCA 
CGT TTCA 

TCT CAT CCG AAG ACG ACA CGT ACC 


iVd = 91 


101 

100* 

99* 

98* 

97* 

96* 

95* 


CAG AAT CCG TGC TAG TAC CTC AAT ATA GAC TCC CTG TAG ATG TTA ATA CTT CTC CTA TGA CTT CTA ATT GAC CCA TGA CTA TCC TCA CCT CCA CCG TTT CA 

AGA ATC CGT GCT AGT ACC TCA ATA TAG ACT CCC TGT AGA TGT TAA TAC TTC TCC TAT GAC TTC TAA TTG ACC CAT GAC TAT CCT CAC CTC CAC CGT TTC A 

GAA TCC GTG CTA GTA CCT CAA TAT AGA CTC CCT GTA GAT GTT AAT ACT TCT CCT ATG ACT TCT AAT TGA CCC ATG ACT ATC CTC ACC TCC ACC GTT TCA 

AAT CCG TGC TAG TAC CTC AAT ATA GAC TCC CTG TAG ATG TTA ATA CTT CTC CTA TGA CTT CTA ATT GAC CCA TGA CTA TCC TCA CCT CCA CCG TTT CA 

ATC CGT GCT AGT ACC TCA ATA TAG ACT CCC TGT AGA TGT TAA TAC TTC TCC TAT GAC TTC TAA TTG ACC CAT GAC TAT CCT CAC CTC CAC CGT TTC A 

TCC GTG CTA GTA CCT CAA TAT AGA CTC CCT GTA GAT GTT AAT ACT TCT CCT ATG ACT TCT AAT TGA CCC ATG ACT ATC CTC ACC TCC ACC GTT TCA 

CCG TGC TAG TAC CTC AAT ATA GAC TCC CTG TAG ATG TTA ATA CTT CTC CTA TGA CTT CTA ATT GAC CCA TGA CTA TCC TCA CCT CCA CCG TTT CA 


Variable-sequence 


73 

TA 

73 

E8A10 

73 

E8A17 

73 

E8A26 

73 

E8A38 

Structural defects 

69“ 

C:C mismatch 

97"" 

Nick 1 

97"" 

Nick 2 

97"" 

Double Nick 


CAG AAT CCG TAG CTC TAG CAC CGC TTA AAC GCA CGT ACG CGC TGT CTA CCG CGT TTT AAC CGC CAA TAG GAT T 
CAG AAT CCG TTT TTA TTT ATC GCC TCC ACG GTG CTG TTT TTT TTT TCT GTT GGC CGT GTT ATC TCG AGT TAG T 
CAG AAT CCG TTT TTA TTT ATC GCC TCC ACG GTG CTT TTT TTT TTT TTT TTT GGC CGT GTT ATC TCG AGT TAG T 
CAG AAT CCG TTT TTA TTT ATC GCC TCC TTT TTT TTT TTT TTT TTT TTT TTT TTC CGT GTT ATC TCG AGT TAG T 
CAG AAT CCG TTT TTA TTT ATC GTT TTT TTT TTT TTT TTT TTT TTT TTT TTT TTT TTT TTT ATC TCG AGT TAG T 


CAG AAT CCG TGC TAG TAC CTC AAT ATA GAC TCC CTT TGA CCC ATG ACT ATC CTC ACC TCC ACC GTT TCA 

CAG AAT CCG TGC TAG TAC CTC AAT ATA GAC TCC CTA I TGT TAA TAC TTC TCC TAT GAC TTC TAA TTG ACC CAT GAC TAT CCT CAC CTC CAC CGT TTC A 

CGA TCA TGG AGT TAT ATC TGA GGG ATA CAA TTA TGA AGA GGA TAC TGA AGA TTA ACT GGG TAC TGA TAG GAG TGG AGG TGG CAA AGT GTC TTA GGC A 

CAG AAT CCG TGC TAG TAC CTC AAT ATA GAC TCC CTA TGT TAA TAC TTC TCC TAT GAC TTC TAA TTG ACC CAT GAC TAT CCT CAC CTC CAC CGT TTC A 

CGA TCA TGG AGT TAT ATC TGA GGG ATA CAA TTA TGA AGA GGA TAC TGA AGA TTA ACT GGG TAC I TGA TAG GAG TGG AGG TGG CAA AGT GTC TTA GGC A 

CAG AAT CCG TGC TAG TAC CTC AAT ATA GAC TCC CTA I TGT TAA TAC TTC TCC TAT GAC TTC TAA TTG ACC CAT GAC TAT CCT CAC CTC CAC CGT TTC A 

CGA TCA TGG AGT TAT ATC TGA GGG ATA CAA TTA TGA AGA GGA TAC TGA AGA TTA ACT GGG TAC I TGA TAG GAG TGG AGG TGG CAA AGT GTC TTA GGC A 


Table SI. Unless otherwise stated,^ cyclization sequences are from WH [S3], with sticky end (Ns) highlighted in red, and 
the remainder of the sequence (N^) in black. Sequences are written 5' to 3^ For conciseness, other than the initial example of 
Nbp = 67, a second strand is omitted where it is the reverse complement of the first strand. 

^ Non-WH sequence. 

™ C:C Mismatch highlighted (blue). 

^ Nicks highlighted with bar | (blue). 


For the computation of jeq and idyn^ we use average 
values for 336 nM. For sequence-dependent results at 
Nbp = 73, a length where the dimerization equilibrium 
constant was computed, we use the appropriate length- 
specific value instead of the average. 


B. Initial stress in cyclized system 

There is a general trend towards a larger activation 
free-energy barrier to cyclization (AGjy^,) at shorter Nbp. 
This is due to the bending stress imposed upon the system 
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Vi/bp 

[Conc]/nM 


Average-base 



Sequence-dependent 




iiTeq” / 10^^ M-1 

/ ksT / knT 

/ 10^^ M 

/ knT AGl^,^ / knT 

20;20 

336 

l.lliO.ll 

14.42 zb 0.11 

27.09 zb 0.10 

0.42 zb 0.05 

14.13 zb 0.07 

25.75 zb 0.11 

57;57 

336 

0.87 zb 0.05 

14.62 zb 0.04 

27.05 zb 0.05 

0.54 zb 0.03 

14.01 zb 0.01 

25.79 zb 0.06 

63;63 

336 

0.78 zb 0.04 

14.69 zb 0.15 

27.00 zb 0.05 

0.33 zb 0.05 

14.19 zb 0.05 

25.60 zb 0.13 

91;91 

336 

0.91 zb 0.15 

14.66 zb 0.07 

27.13 zb 0.16 

0.44 zb 0.03 

14.20 zb 0.04 

25.86 zb 0.08 

Avg 

336 

0.92 zb 0.06 

14.60 zb 0.06 

27.07 zb 0.05 

0.43 zb 0.03 

14.12 zb 0.03 

25.75 zb 0.05 

10;10 

2690 

1.22 zb 0.19 

12.54 zb 0.05 

27.38 zb 0.16 

0.46 zb 0.02 

12.11 zb 0.11 

25.91 zb 0.04 

28;29 

2690 

1.11 zb 0.09 

12.66 zb 0.10 

27.41 zb 0.08 

0.82 zb 0.12 

12.60 zb 0.04 

27.03 zb 0.15 

30;33 

2690 

1.11 zb 0.17 

12.64 zb 0.12 

27.40 zb 0.15 

1.15 zb 0.12 

12.53 zb 0.05 

27.32 zb 0.11 

45;46 

2690 

1.04 zb 0.09 

12.55 zb 0.08 

27.24 zb 0.09 

0.48 zb 0.05 

12.08 zb 0.04 

25.90 zb 0.11 

Avg 

2690 

1.12 zb 0.07 

12.60 zb 0.05 

27.36 zb 0.06 

0.73 zb 0.08 

12.32 zb 0.06 

26.54 zb 0.15 


Table S2. Equilibrium constant for dimerization and activation free-energy barriers to dimerization and undi¬ 
merization as defined in the main text for both the oxDNA average-base and sequence-dependent parameterizations. 

Simulations were performed at T = 298 K. The first column gives the lengths of duplex regi ons of the two monomers. The 
complementary sticky ends are of identical sequence and length Ns = 10. Sequences from Table (for [Cone] =336 nM the rel¬ 
evant sequences are those with iV^p = iV^ + lO; for [Cone]=2690 nM they are those with iV^p equal to the total length of the dimer). 


by the formation of the initial base pair (Qpp = !)• 
very shortest lengths 30 — 45 , this trend is reversed 

because the complementary single-stranded sticky ends 
(of length Ns = 10) are sufficiently long relative to TVq 
that an initial base pair can form with less bending stress 
in the duplex region. States at the top of the free-energy 
barrier (i.e. Qbp = 1) for cyclization and dimerization 
at TVbp = 30,101 are compared in Figure SI, clearly 
showing bending in the cyclized molecules, but not the 
dimerized molecules. Importantly, at TV^p = 30, the 
relatively long single-stranded region {Ns = ^/2 x Vd = 10) 
reduces the requirement to bend the duplex (FigurejS^c)), 
compared to the strong bending necessary at A/^p = 101 
(Figure [sTta)). 


C. Comparison to FRET experiments 


1. Length of complementary sticky ends 


We compare oxDNA to available V&H equilibrium data 
[S3] for Aq = 91, As = 8 — 10. Specifically, V&H report 
the fractional occupancy of the cyclized state /eye- Akin to 
the nomenclature in Section |S1 C| /eye +/open = 1, where 
/open is the fractional occupancy of the open (uncyclized) 
state. 

Fractional occupancy can be a somewhat misleading 
measure by which to compare systems, because the en¬ 
tire transition region /eye ~ 0.25 — 0.75 is covered by a 
free-energy difference of ^ 2 while a similar free- 

energy difference at the extrema (/eye —> 0 or 1) would 
be indistinguishable. We therefore make our equilibrium 


comparison in terms of AGcyc, where 


AGcyc = -A/B^In 

(S8) 

and 


T^cyc /eye /eye 

^eq - f ~ 1- f ' 

Jopen -L Jcyc 

(S9) 


OxDNA appears to be in good agreement with experi¬ 
ment; although, this may be partly coincidental as V&H 
do not report either the temperature (presumed to be 
298 K) or a salt concentration (V&H give their imaging 
buffer cation concentration as [Na^j = 500 — 1000 mM or 
[Mg^^j = 10—30 mM; recall that oxDNA is parameterized 
to [Na^j = 500 mM). At Ag = 9 and Ag = 10, oxDNA 
and experiment differ by ^ 1 T, while at Ag = 8, 
oxDNA under-reports experiment by ^ 3 ^ (cyclized 

state is less likely in oxDNA than experiment). 

As expected from the SantaLucia model [S3 ED of 
nearest-neighbour thermodynamics, the linear relation¬ 
ship between oxDNA results in Figure [^(a) is reffect- 
ive of the ^ 3 A;p T stabilization of each additional base 
pair. A possible explanation for the discrepancy between 
oxDNA and V&H |S3| is that hairpin formation in the Ag 
region could decrease the fraction of cyclized molecules; 
but given the sequence, hairpin formation should be neg¬ 
ligible. In particular, we cannot explain the discrepancy 
between oxDNA and the V&H data for Ag = 8, although 
we note that identifying a very low yield in experiment 
can be challenging, for instance due to the presence of 
impurities inducing large changes in AGcyc- 
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Cyclized 
(a) JVbp = 101 



(c) iVbp = 30 


Dimerized 

(b) 



(d) 




Figure SI. OxDNA representations of Qpp = 1 states for cyc¬ 
lized and dimerized configurations, highlighting the bending 
required to form the first base pair. Dimerized configurations 
in (b) and (d) are shown for iV^i + iVd 2 + iVg = instead 
of length A/j^p monomers. For (b), iV^i = 45, N(^2 — 46 and 
Ns = 10. For (d), N^i = N^2 = 10 and Ns = 10. 


2. Nicks and mismatches 


We examine the role of structural defects, namely nicks 
and mismatches, on cyclization using the average-base 
oxDNA parameterization. Unsurprisingly, we find that 
nicked and mismatches molecules cyclized more readily 
than their intact counterparts (Figure [S^. 

V&H briefly examine the impact of nicks and mis¬ 
matches on the rate of cyclization, recording /eye versus 
time (Figure S3 in Ref. [S3]). It does not appear that 
all molecules have reached an equilibrium plateau in /eye, 
therefore a concrete equilibrium comparison to oxDNA 
is difficult. We do, however, observe good agreement 
between the oxDNA and the experimental observations 
of Vfc H re garding the stabilizing effect of various motifs 
(Table [S^. However, we do not reproduce the unexpec¬ 
ted similarity between the experimental values for intact 
sequences at A^p = 69 and A^p = 97. This is consistent 
with our inability to reproduce V&H’s j-factors for their 
shortest sequences. 


(a) 


Fh 


<] 
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10 
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-5 
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• oxDNA Avg-seq - 

• oxDNA Seq-dep 
■ V&H 
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Qbp 


Figure S2. Dependence of cyclization thermodynamics on Ag. 

(a) Comparison between oxDNA and V&H [S3] for AGcyc, 
the free energy of the cyclized compared to the open state. 

(b) Free-energy profiles for cyclization. The reference state 
corresponding to AC = 0 is set to be at Qpp = 1. All 
results are for A^j = 91, and using the oxDNA average-base 
parameterization, unless stated. 



OxDNA 

* 

Experiment 


T^cyc 

-^eq 

AAG/keT^ 

AAG/keT^ 

iVbp = 69 

Intact 

0.029 ± 0.007 


1.10 

Mismatch 

0.052 ±0.016 

-0.6 

3.14 -1.1 

II 

Intact 

2.99 ±0.35 


0.93 

Nick 1 

18.5 ±0.9 

-1.8 

1.73 -0.6 

Nick 2 

16.2 ± 1.0 

-1.7 


Double Nick 

50.3 ±2.4 

-2.8 

8.78 -2.3 


Table S3. for oxDNA and experiment (V&H |S3|b 

The experimental assays were not intended to probe the equi¬ 
librium behaviour of the systems; as such, the last point 
of kinetic cyclization versus time assay may not represent 
equilibrium values. To emphasize the qualitative difference 
between intact molecules and those with structural defects 
(nicks, mismatches), oxDNA results are for the average-base 
parameterization. 

Converted from yield /eye to using Eq. (S9). 

^ A AG is computed as the difference in the AC of cyclization 
between intact states and those with structural defects (nick 
or mismatch). 
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Qbp 

Figure S3. Free-energy profiles for cyclization for (a) an 
intact duplex compared to a duplex with a single mismatch at 
iVbp = 69, and (b) intact versus nicked duplexes at iV^p = 97, 
using the oxDNA average-base parameterization and Ns — 10 
for all sequences (Table [ST| . 


D. Comparison to ligase experiments 

Comparison to the ligase experiments of Du et ai [S9] 
and Cloutier & Widom (C&W) |S10| highlight some in¬ 
teresting features of the oxDNA model (Figure [S^. For 
example, while both the Du et al. experimental and 
oxDNA model results are similar to the SY WLC ex¬ 
pression for TVbp > 100, there is a signihcant absolute 
deviation in j-factor due to differences in the relevant 
mechanical properties for oxDNA and those seemingly 
appropriate for the conditions of the experiment of Du et 
al. Specihcahy, they use NEB T4-hgase buffer, which has 
a salt concentration of [Mg^^] = 10 mM, while oxDNA 
is parametrized to [Na^] = 500 mM. Du et a/.’s ht to 
their data gives a weaker torsional stiffness (2.4 x 10~^^ 
versus 4.75 x 10~^^ Jm), longer persistence length (47 
versus 41.82 nm) and longer pitch length (10.54 versus 
10.36 bp/turn) than for oxDNA. The comparison also 
highlights how relatively small deviations in DNA mech¬ 
anical properties (in particular a ^ 10 % change in the 
persistence length) may shift the apparent jeq by an order 
of magnitude. The differences in SY WLC expression 
parameters are reasonable given the different solution 
conditions (the persistence length decreases with increas¬ 
ing salt concentration [Sill IS12| and temperatures (e.g. 
jeq is enhanced by a factor of ~ 10 between 5 °C and 
42 °C |S13| b We note that the oxDNA persistence length 
is reasonable for the high salt conditions used by V&H 



Ybp 

Figure S4. OxDNA comparison to the experiments of Du et 
al. [S9] and Cloutier & Widom [S10| . (a) Du et al. equilibrium 
j/’-factor via a ligase-based assay (green circles) compared 

to the SY WLC expression with Du et al. parameters 

(gray dotted) and the oxDNA model (gray marks). 

The top envelop of with oxDNA parameterized is dashed 

(gray dashed) (b) C&W ligase-based j-factor compared to 
Myn^^ (black marks) and^°q^^^ (gray marks). The maxima 
envelope of ^ is given with a gray dashed line. 

EE]. 

Du et al. only investigated lengths near the maxima in 
ieq^^' ^ consequence, their experiments are unable to 
conhrm our prediction of a decrease in the magnitude of 
the oscillations in which we attribute to the af¬ 

fect of teardrop conhgurations. Whether the experiments 
could potentially see the effects of teardrop conhgurations 
also depends on how the ligase activity depends on the 
nature of the DNA conhrmation. Indeed, Vologodskii et 
al. mention the possibility of poor ligation efficiency for 
teardrop conhgurations m\- We also note that because 
of the potential role of teardrop conhgurations, the tor¬ 
sional modulus extracted from hts to the SY expression 
should be interpreted with caution, and it is interesting 
to note that the values typically obtained from these hts 
[SbllS^IST^ are signihcantly lower than those from single¬ 
molecule experiments with magnetic tweezers |S16[ IS17] . 

Since it is likely that C&W conducted their experiment 
at too high a ligase concentration, it is difficult to char¬ 
acterize their j-factor in terms of either or jeq. To 
isolate high ligase concentration as the culprit behind an 
apparent deviation from Du et al. used conditions 

identical to C&W, namely buffer ([Mg^^] = 10 mM) and 
temperature (T = 21 °C for Du et al.^ T = 20°C for 
C&W and T = 25°C for oxDNA). Intriguingly, C&W’s 
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Figure S5. The ensemble average of the oxDNA capture 
radius (i?capture) as a function of length iV^p at fixed Ns = 10. 
Error bars represent the standard deviation. 


results mainly lie inbetween 


E. Capture radius 


Comparing results for with a WLC model requires 
that the finite dimensions of the sticky ends be taken into 
account. Defining a “capture radius”, a measure of the 
distance within which the sticky ends can hybridize, is 
a common way of accounting for the finite size of DNA 
sticky ends. As oxDNA naturally accounts for the finite- 
size of the sticky ends, we have the ability to test the 


capture radius assumption. 

In oxDNA, we define the capture radius i^capture as the 
distance between the first and last bases of the duplex 
region at = 1 (Figure Sl(a,c)). Nucleotide 

positions are used to report distance, instead of the mid¬ 
point along the helical axis, because fraying may occur at 
the ends of the duplex. The choice of nucleotide strand 
does not impact our results. 

Our ensemble average capture radius, (^capture) is 
shown in Figure (^capture) is roughly constant at 
^bp ^ 100 — 200, increases slowly between A^p « 
50 — 100, and then more rapidly for A^p « 40 — 50, 
before decreasing for A^p < 37. The saturation value 
((^capture) ~ 4nm) is quite close to the capture radius of 
5nm that has been used to interpret V&H’s experiments 
with 10-base sticky ends |S3[ IS18| . 

The free-energy cost associated with forming the ini¬ 
tial contact is partitioned between that for bending the 
duplex and that for stretching the single-stranded sticky- 
ends. At shorter A^p, the amount of bending required 
for initial binding increases, as does the stretching of the 
single-stranded sticky-ends. At short lengths, this trend 
is reversed when the length of the unperturbed duplex 
approaches the capture radius at = 1. Interestingly, 
the distribution of ^capture is roughly Gaussian, with 
similar standard deviation for all lengths. 

While we explicitly disallow misbonding, allowing mis- 
bonding may slightly increase the capture radius, an ef¬ 
fect that should be more prominent for longer sticky-ends. 
Overall, the 5 nm value for the capture radius used in the 
interpretation of V&H’s experiments appears not unreas¬ 
onable. 
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